Load packages
library(Seurat)
library(ggplot2)
library(cowplot)
Load previously saved datasets and add annotation
cdS1 <- readRDS("~/Documents/Experiments1/scRNAseq/Cerebellum_publication/allCells/data/Seurat3.rds")
cdS1@meta.data$dataset <- "E13"
colnames(cdS1@meta.data)[1:2] <- c("nCount_RNA","nFeature_RNA")
cdS2 <- readRDS("./E14/Seurat_E14.rds")
cdS3 <- readRDS("./E15/Seurat_E15.rds")
cdS4 <- readRDS("./E16/Seurat_E16.rds")
cdS1@meta.data$dataset1 <- "E13"
cdS2@meta.data$dataset1 <- "E14"
cdS3@meta.data$dataset1 <- "E15"
cdS4@meta.data$dataset1 <- "E16"
Dataset preprocessing
cdS <- merge(cdS1, c(cdS2, cdS3, cdS4))
cb.list <- SplitObject(object = cdS, split.by = "dataset1")
Normalize data and find variable genes
for (i in 1:length(x = cb.list)) {
cb.list[[i]] <- NormalizeData(object = cb.list[[i]], verbose = FALSE)
cb.list[[i]] <- FindVariableFeatures(object = cb.list[[i]],
selection.method = "vst", nfeatures = 2000, verbose = FALSE)
}
Integrate the deatasets
anchors <- FindIntegrationAnchors(object.list = cb.list, dims = 1:30)
cdS.i <- IntegrateData(anchorset = anchors, dims = 1:30)
switch to integrated assay
DefaultAssay(object = cdS.i) <- "integrated"
Run the standard workflow for visualization and clustering
cdS.i <- ScaleData(object = cdS.i, verbose = FALSE)
cdS.i <- RunPCA(object = cdS.i, npcs = 40, verbose = FALSE)
ElbowPlot(cdS.i, ndims = 40, reduction = "pca")
Find clusters
cdS.i <- FindNeighbors(cdS.i, reduction = "pca", assay = "integrated", dims = 1:23)
Computing nearest neighbor graph
Computing SNN
cdS.i <- FindClusters(cdS.i, modularity.fxn = 1,
resolution = 1, algorithm = 1, n.start = 10, n.iter = 10,
random.seed = 0, temp.file.location = NULL, edge.file.name = NULL,
verbose = F)
Inspect cell clusters
myplots <- list()
for (n in 1:length(res)){
res1 <- paste("integrated_snn_res",res[n],sep=".")
Idents(object = cdS.i) <- cdS.i@meta.data[,res1]
myplots[[n]] <- DimPlot(cdS.i, reduction = "tsne", group.by = "ident", label=T)+
NoAxes()+NoLegend()
}
myplots1 <- list()
for (n in 1:length(res)){
res1 <- paste("integrated_snn_res",res[n],sep=".")
Idents(object = cdS.i) <- cdS.i@meta.data[,res1]
myplots1[[n]] <- DimPlot(cdS.i, reduction = "umap", group.by = "ident", label=T)+
NoAxes()+NoLegend()
}
gridExtra::grid.arrange(grobs= myplots, ncol = 2, top = paste0("resolutions: ", res))

gridExtra::grid.arrange(grobs= myplots1, ncol = 2, top = paste0("resolutions: ", res))

res1 <- paste("integrated_snn_res",2,sep=".")
Idents(object = cdS.i) <- cdS.i@meta.data[,res1]
Visualize the results with t-SNE

Visualize the results with UMAP

Visualize the results with UMAP
cdS.i <- RunUMAP(object = cdS.i, dims = 1:23, reduction = "pca",
features = NULL, assay = "integrated", nneighbors = 65L, max.dim = 2L,
min.dist = 0.4, reduction.name = "umap", reduction.key = "UMAP_",
metric = "correlation", seed.use = 42)
DimPlot(cdS.i, reduction = "umap", group.by = "ident", label=T)
DimPlot(cdS.i, reduction = "umap", group.by = "cellType", label=T, split.by = "dataset1")+
NoAxes()+NoLegend()

DimPlot(cdS.i, reduction = "tsne", group.by = "cellType", label=T, split.by = "dataset1")+
NoAxes()+NoLegend()
